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Abstract: This paper presents a sensor system for robot localization based on the information 
obtained from a single camera attached in a fixed place external to the robot. Our approach 
firstly obtains the 3D geometrical model of the robot based on the projection of its natural 
appearance in the camera while the robot performs an initialization trajectory. This paper 
proposes a structure-from-motion solution that uses the odometry sensors inside the robot as 
a metric reference. Secondly, an online localization method based on a sequential Bayesian 
inference is proposed, which uses the geometrical model of the robot as a link between image 
measurements and pose estimation. The online approach is resistant to hard occlusions and 
the experimental setup proposed in this paper shows its effectiveness in real situations. The 
proposed approach has many applications in both the industrial and service robot fields. 

Keywords: robotics; vision sensor; intelligent spaces 



1. Introduction 

This paper presents a sensor system composed of a single camera attached to a fixed position and 
orientation in a bounded environment ( indoor workplace) which is observing a mobile robot. The aim 
of the sensor system is to obtain the orientation and position (i.e., pose) of a mobile robot using both 
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visual information retrieved by the camera and relative odometry readings obtained from the internal 
sensors of the robot. The camera acquisition and image processing tasks are executed in a specialized 
hardware, which also controls the behavior and internal sensors of the mobile robot through a wireless 
channel. The proposed schema allows the robot to perform complex tasks without requiring dedicated 
processing hardware on it. This approach is sustained in the Intelligent Space [1] concept and it can 
be equally applied to multiple scenarios, specially both in the industrial field (e.g., automatic crane 
positioning, autonomous car parking) and in the service fields (e.g., wheelchair positioning in medical 
environments or autonomous driving of mobile platforms ). The single camera solution presented in 
this paper allows to cover large areas with less cameras compared to multiple camera configurations 
where overlapped areas are mandatory. This feature reduces the cost and improves the reliability of the 
intelligent space philosophy. 

In this paper, we suppose that the camera is correctly calibrated and thus the parameters governing the 
projection model of the camera are previously known. To connect the pose of the robot with information 
found in the image plane of the camera, we propose to obtain a 3D geometric model of the mobile robot. 
Such model is composed of several sparse points whose coordinates represent some well-localized points 
belonging to the physical structure of the robot. These points are determined by image measurements, 
called image features, which usually correspond to corner-like points due to texture changes or geometry 
changes such as 3D vertexes. Usually in industrial fields the image features are obtained by including 
some kind of artificial marker on the structure of the robot ( infrared markers or color bands). These 
methods are very robust and can be used to recognize human activity and models with high degrees of 
freedom (AICON 3D online [2] or ViconPeak online systems [3]). However, in this paper we want to 
minimize the required "a priori" knowledge of the robot, so that it is not necessary to place artificial 
markers or beacons on it to detect its structure in the images. 

Independent of the nature of the image features, the information obtained from a camera is naturally 
ambiguous and thus some sort of extra metric information are required in order to solve such ambiguity. 
In this work, we propose to use the odometry sensors inside the robot {i.e., wheel velocities) to act as the 
metric reference. 

The general diagram of the algorithm proposed in this paper is shown in Figure 1 . It shows a clear 
division of the processes involved in obtaining the pose of the robot: first we denote as "Initialization of 
Pose and Geometry" to those processes necessary to start up the system, such as the 3D model of the 
robot and the initial pose it occupies. The initialization consists of a batch processing algorithm where 
the robot is commanded to follow a certain trajectory so that the camera is able to track some points of 
the robot's structure under different viewpoints jointly with the recording of the odometry information. 
All this information is combined to give the 3D model of the robot and the initial pose it occupies. 

Given the initialization information, the second group of processes, named "Sequential Localization", 
provides the pose of the robot in a sequential manner. It is composed of a pose estimator, given 
odometry readings and a pose correction block which combines the estimation of the pose with image 
measurements to accurately give a coherent pose with the measurements. This algorithm operates 
entirely on-line and thus the pose is available at each time sample. 

Both group of processes are supplied with two main sources of information: 

1. Image measurements: they consist of the projection in the camera's image plane of certain 
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points of the robot's 3D structure. The measurement process is in charge of searching coherent 
correspondences through images with different perspective changes due to the movement of 
the robot. 

2. Motion estimation of the robot: The odometry sensors built on-board the robot supply the 
localization system with an accurate motion estimation in short trajectories but that is prone to 
accumulative errors in large ones. 



Figure 1. General diagram of the proposed localization system using a vision sensor and 
odometry readings. 
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1.1. Previous Works 

Despite the inherent potential of using external cameras to localize robots, there are relative few 
attempts to solve it compared to the approach that consider the camera on-board the robot [4, 5]. 
However, some examples of robot localization with cameras can be found in the literature, where the 
robot is equipped with artificial landmarks, either active [6, 7] or passive ones [8, 9]. In other works 
a model of the robot, either geometrical or of appearance [10, 11], is learnt previously to the tracking 
task. In [12, 13], the position of static and dynamic objects is obtained by multiple camera fusion inside 
an occupancy grid. An appearance model is used afterwards to ascertain which object is each robot. 
Despite the technique used for tracking, the common point of many of the proposals found in the topic 
comes from the fact that rich knowledge is obtained previously to the tracking, in a supervised task. 

In this paper the localization of the robot is obtained in terms of its natural appearance. We propose a 
full metric and accurate system based on identifying natural features belonging to the geometric structure 
of the robot. On most natural objects we can find points whose image projection can be tracked in the 
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image plane independently of the position the object occupies and based on local properties found in the 
image (i.e., lines, corners or color blobs). Those points are considered natural markers, as they serve as 
reference points in the image plane that can be easily relate with their three-dimensional counterparts. 
The set of methods focused on tracking natural markers have become a very successful and deeply 
studied topic in the literature [14, 15], as they represent the basic measurements of most of existing 
reconstruction methods. 

Scene reconstruction from image measurements is a classic and mature discipline in the computer 
vision field. Among the wide amount of proposals it can be highlighted those grouped under the name 
"Bundle Adjustment" [16, 17]. Their aim is essentially to estimate, jointly and optimally, the 3D 
structure of a scene and the camera parameters from a set of images taken under some kind of motion. 
(i.e., it can be the camera that moves or equally some part of the scene w.r.t. the camera). 

In general terms, Bundle Adjustment reconstruction methods are based in iterative optimization 
methods which try to minimize the image discrepancy between the measured positions of the 3D model 
and the expected ones using the last iteration solution. The discrepancy balances the contribution of the 
measurements into the final solution and plays an important role in this paper. Our main contribution 
is a redefinition of the discrepancy function using a Maximum Likelihood approach which takes into 
account the statistical distribution of the error. This distribution is especially affected by the odometry 
errors which are accumulative in long trajectories. 

On the other hand, once a geometric model is obtained using a structure-from-motion approach, 
its pose with respect to a global coordinate origin can be easily retrieved by measuring the projection 
of the model in the image plane. This problem, commonly known as the Perspective n Point Problem 
(PnP), has received considerable attention in the literature, where some accurate solutions are found such 
as [18] or the recent global solution proposed in [19]. In this paper we instead follow a filtering approach, 
where not only image measurements but also last time pose information and odometry information are 
used to obtain the pose. This approach, which is based on the use of a Kalman Filter, is much more 
regular than solving the PnP problem at each time instant and will be described in this paper. 

The paper is organized as follows. In Section 2. the notation and mathematical elements used in the 
paper are described. In Section 3. the description of the initialization algorithms is given. The online 
Kalman loop is explained in Section 4. and some results in a real environment are shown in Section 5. 
Conclusions are in Section 6. 

2. Definitions and Notation 

This section presents a description of the symbols and mathematical models used in the rest of 
the paper. 

2. 1. Robot and Image Measurements 

The robot's pose at time k is described by a vector X^. We suppose that the robot's motion lie on 
the plane 2 = 0 referred from the world coordinate origin Ow (See Figure 2). The pose vector X k 
is described by 3 components (xk,yk,&k), corresponding to two position coordinates Xk,Vk and the 
rotation angle oi^ in the z axis. The motion model Xk = g(Xk-i, Uk) defines the relationship between 
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the current pose X k with respect to the previous time one X k _i and the input U k given by odometry 
sensors inside the robot. In our proposal the robot used is a differential wheeled platform and thus 
U k = (Vl k ,Vt k ) T , where Vl k and Q k describe the linear and angular speed of the robot's center of 
rotation (Or in Figure 2) respectively. The motion model g is then very simple in function of the 
discretized linear speed Vl k and the angular speed fl, k . 

The robot's geometry is composed of a sparse set of N 3D points M. = {M 1 , • • • , M N } referred from 
the local coordinate origin Or described by robot's pose Xk. The points M l are static in time due to 
robot's rigidness, and thus, no temporal subindex is required for them. Function M l x = t(Xk, M l ) 
uses actual pose X k to express M % in the global coordinate origin Ow that X k is referred to (see 
Figure 2): 



M Xk =t(X k ,M t )=R k M* + T k 



(1) 



where: 
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Figure 2. Spatial relationship between world's coordinate origin Ow, robot's coordinate 
origin Or and camera's coordinate origin Oq. 
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The augmented vector X%, which is the state vector of the system, is defined as the concatenation in 
one column vector of both the pose X k and the set of static structure points M. : 



X a k = ({X k ) T (MY ... (M^f 
An augmented motion model g a is defined as the transition of the state vector: 
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It must be noticed that the motion model g a leaves the structure points contained in the state vector 
untouched, as we suppose that the object is rigid. 

The whole set of measurements obtained in the image plane of the camera at time k is denoted by Y k . 

Each measurement inside Y k is defined by a two-dimensional vector y k = (u l k v^) , describing the 
projection of each point inside M. in the image plane, where i stands for the correspondence with the 
point M l . It must be remarked that, in general, Y k does not contain the image projection of every point 
in the set M as some of them can be occluded depending on the situation. 

The camera is modeled with a zero-skew "pin-hole" model (see [17] for more details), with the 
following matrix K c containing the intrinsic parameters: 



(5) 



The extrinsic parameters of the camera (i.e., position and orientation of the camera w.r.t. Ow) are 
described using the rigid transformation R c , T c (i.e., rotation matrix and translation vector). The matrix 
R c is defined by three Euler angles a c , (3 C , %. The vector T c can be decomposed into its three spatial 
coordinates T x ,T y ,T z . All calibration parameters can be grouped inside vector P: 

P = (fu, fv, u 0 , v 0 , a c , p c , 7c, T x , T y , T Z ) T (6) 

Given a single measurement y l k from the set Y k , its 2D coordinates can be expressed using the 
"pin-hole" model and the aforementioned calibration parameters: 
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\\K C (R c M l Xk + T c ) = \\K C (R c R k M* + R c T k + T c ) (7) 



where A^. represents a projective scale factor that can be obtained so that the third element of the right 
part of Equation (7) is equal to one. It is important to remark that the projection model, although simple, 
depends in a non-linear fashion w.r.t. M % due to the factor A^,. For compactness the projection model of 
Equation (7) can be described as the function h: 

yl = h(X k ,M\P) (8) 
In the same way, the whole vector Y k can be expressed with the following function: 

Y k = h a (X a k ,P)=(h(X k ,M\P) T ■■■ h(X k ,M N ,P) T Y (9) 
2.2. Random Processes 

This paper explicitly deals with uncertainties by describing every process as a random variable with 
a Gaussian distribution whose covariance matrix stands for its uncertainty. The random processes are 
expressed in bold typography (i.e., Xk is the random process describing the pose and X k is a realization 
of it) and each of them are defined by a mean vector and a covariance matrix. Therefore is defined 
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by its mean X% and covariance (or simplifying J\f(X%, ££)). The following processes are considered 
in the paper: 

• PoseX k = Af(X k , Efc) and 3D model M = J\f(M, Em) processes. Its joint distribution is encoded 
inXJ=^(X^Eg). 

• Measurement process Yk = ftf(Y~k, EyJ, whose uncertainty comes from errors in image detection. 

• Odometry input values Uk = N(Uk, £<i/ fc )> which are polluted by random deviations. 
3. Initialization Process 

The initialization step consists of obtaining the initial pose the robot occupies and a 3D model of 
its structure using visual features (i.e., to obtain an initial guess of Xq). The importance of this step 
is crucial, as the robot's geometric model serves as the necessary link between robot's pose and image 
measurements. The initialization allows to use afterwards the online approach presented in Section 4. 

A delayed initialization approach is proposed, based on collecting image measurements and the 
corresponding odometry position estimation along a sequence of time (i.e., k — 1, • • • , K). After taking 
a sufficient amount of measurements of the object in motion, an iterative optimization method is used to 
obtain the best solution for Xq according to a cost criterion. The odometry readings are used in this step 
as metric information, which allows to remove the natural ambiguity produced by measurements from 
a single camera. The main problem of using odometry as a metric reference for reconstruction comes 
from the accumulative error it presents, which is a very well-known problem in dead-reckoning tasks. 



Figure 3. Maximum likelihood initialization by means of reducing the residual of the 
expected measurements (red diamonds) and the measured image trajectories (blue circles). 
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The initialization algorithm consists of a Maximum Likelihood (M.L.) estimation of the vector 
Xq, which does not estimate the initialization value itself but its equivalent Gaussian distribution 
= A/'(X 0 1 ,Eg) by using a "Bundle Adjustment" method (See Figure 3). The M.L. approach 
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(see [16, 17] for more details) allows to properly tackle the growing drift present in the odometry 
estimation by giving fair less importance or weight to those instants where the odometry is expected 
to be more uncertain. 

In the following sections, the initialization algorithm is explained in detail. 

3.1. Odometry Estimation 

A set of odometry readings, Ui, - ■ ■ , Uk are collected over a set of K consecutive time samples 
k = 1, • • • , K, which corresponds to the initialization trajectory performed by the robot. Using a motion 
model g given initial position X 0 , an estimation of the robot's pose in the whole initialization sequence 
is obtained as follows: 

Xi = ff(X 0 ,Ui) 
X 2 = ff (X 1; U 2 ) 

X K = <?(X K -l,Uk) (10) 

where we recall that Xk and Uk denote Gaussian processes. Using expression (10), and by propagating 
the statistical processes through function g, the joint distribution p(Xi, ■ ■ ■ , X K \X 0 ) can be obtained. 
The propagation of statistical processes through non-linear functions is in general a very complex task, 
as it requires to solve integral equations. However in this paper we will follow a first order approximation 
of the non-linear functions centered in the mean of the process (see [20] for more details). Using this 
technique the Gaussian processes Xk and Uk can be iteratively propagated through function g at the cost 
of being an approximation. This paper shows that despite being an approximation this approach end in 
a good estimation of the initialization vector. 

By denoting as X to the joint vector containing all poses X = (Xf, ■ ■ ■ , X]^) T , the joint distribution 
of all poses, given the initial state X Q , p(X\X 0 ) is approximated as Gaussian p.d.f. with mean X and 
covariance matrix Ex- 

3.2. Image Measurements 

In this step we describe how to collect the position in the image plane of different points from the 
robot's structure during the initialization trajectory. We base our method in the work [15] where the 
SIFT descriptor is introduced. The process used in the initialization is composed of two steps: 

• Feature -based background subtraction. 

We describe the static background of the scene using a single image, from which a set of features, 
— {vli ' ' ' > Ub Nb } an d its correspondent SIFT descriptors T> 1 = {d\, ■ ■ ■ , df d } are obtained. 
The sets T>b and Tb are our feature-based background model. 

Given an input image, namely I k , we find the sets V k and T k . We consider that a feature 
yl G Fk,d\ G V k belongs to the background if we can find a feature y J b G Tb, d{ G V b such 
that \y l k — y J b \ < R ma x and \d\ — d J b \ < d max . This method although simple shows to be very 
effective and robust in real environments (See Figure 4 ). 
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Figure 4. Feature-based background subtraction method. 




• Supervised tracking algorithm. 

Given the first image of the initialization trajectory, namely I\, we obtain the set of iV features 
T\, once they are cleaned by the aforementioned background subtraction method. We propose a 
method to track them in the whole initialization trajectory. 

By only using SIFT descriptors and its tracking in consecutive frames does not produce stable 
tracks, specially when dealing with highly irregular objects where many of the features are 
due to corners. We thus propose to use a classical feature tracking approach based on the 
Kanade-Lucas-Tomasi (KLT) algorithm [14]. To avoid degeneracy in the tracks, which is a 
very common problem in those kind of trackers, we use the SIFT descriptors to remove those 
segments of the tracks that clearly do not correspond to the feature tracked. This can be done 
easily by checking distance between the descriptors in the track. The threshold used must be 
chosen experimentally so that it does not eliminate useful parts of the tracks. In Figure 5a we 
can see the tracks obtained by the KLT tracker without degeneracy supervision. In Figure 5b the 
automatically removed segments are displayed. 

Figure 5. Comparison between the KLT tracker and the SIFT-supervised KLT version. 




algorithm 

3.3. Likelihood Function 



Using the feature-based algorithm proposed before a set of measurements in the whole trajectory, 
Yi,'" , Y K , is obtained, where each vector Y k contains the projection of iV points from robot's structure 
at time sample k. The set of N points, M 1 , • • • , M N , jointly with the initial pose X 0 , represent the 
initialization unknowns. As the robot's motion does not depend of its structure, the following statistical 
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independence is true: 



p(X ir -- ,X K \X 0 )=p{X li ... ,X K \X 0 ,M\ 
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, X k \Xq) 



(11) 



If all trajectories are packed into vector Y L , the following function put in relation Y L with the 
distribution of expression (11): 



/yA / h a (xi,p) + v 1 \ 



\Y K J \h a (X^P)+V K J 



(12) 



where V k represents the uncertainty in image detection. Using distribution showed in (11), and 
propagating statistics through function (12), the following likelihood distribution is obtained: 



p(Y L \X 0 ,M\-.. ,M N )=p{Y L \XS) 



(13) 



The likelihood function (13) is represented by a Gaussian distribution using a first order approximation 
of (12). It is thus defined by a mean Y L and a covariance matrix E L : 



p(Yl\xs) 



3.4. Maximum Likelihood Approach 



1 



_ e ^Y L -Y L ) T ^- L \Y L -Y L ) 
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(14) 



The likelihood function (14), parametrized by its covariance matrix T, L and its mean Y L , is dependent 
of the conditional parameters and unknown values X§ . 

The "Maximum Likelihood" estimation consists of finding the values for Xq that maximize the 
likelihood function: 



max piYAX^M 1 , 

Xo.Af 1 ,- ,M N 



, M N ) 



(15) 



In its Gaussian approximation and by taking logarithms, it is converted into the following 
minimization problem: 



min {Y L -Y L ) T Y.l\Y L -Y L ) 



(16) 



where Y L is the realization of the process (i.e., the set of measurements from the image) and Y L are 
the expected measurements given a value of the parameters X 0 , M 1 , • • • , M N . 

The configuration of T. L is ruled by the expected uncertainty in the measurements and the statistical 
relation between them. Usually all cross -correlated terms of S L are non-zero, which has an important 
effect in the sparsity of the Hessian used inside the optimization algorithm and consequently in the 
computational complexity of the problem (see [21] for more details). 

The covariance matrix can be approximated assuming either independence between time samples 
(discarding cross-correlation terms between time indexes) or total independence between time and 
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each measurement (discarding all cross-correlation terms except the 2x2 boxes belonging to a single 
measurement). In Table 1, the different cost functions are shown depending on the discarded terms of 
Til. The different approximations of S L have a direct impact in the accuracy and the reconstruction 
error, and the results will be shown in Section 5. 



Table 1. Different cost functions depending on the approximations of T L . 



Type of Approximation 


cost function 


Complete Correlated matrix (M.C.) 


(Y L -Y L fEl\Y L -Y L ) 


2N x 2N block approximation of E L (M.B.) 


Zk=i(Yk-Y k ) T ^(Y k -Y k ) 


2x2 block approximation of S L (M.D.) 


T.tiT^=M-yl) T ^-hyl-yl) 


Identity approximation of E L (M.I.) 


J2k=i J2i=i(yl ~ yk) T (yl ~ vD 



Intuitively, if Sl results to be a identity matrix, the cost function is reduced to a simple image residual 
minimization extensively used in Bundle Adjustment techniques, where in principle all cost differences 
yl — y l k have equal importance: 

N Ki 
0 i=l k=l 

The result of minimizing Equation (16) instead of the usual (17) show significant improvements in 
the reconstruction error. In Section 5. a comparison is shown between the initialization accuracy under 
the different assumptions of the matrix S^, from its diagonal version to the complete correlated form. 
The minimum of (16) is obtained using the Levenberg-Mardquardt [17] iterative optimization method. 

We suppose that during the measurement step there is a low probability of encountering outliers in the 
features. This argument can be very optimistic in some real configurations where more objects appear in 
the scene and produce occlusions or false matchings. For those cases, all cost functions presented in this 
section can be modified to be robust against outliers by using M-Estimators. We refer the reader to [17] 
for more details. 

3.5. Initialization before Optimization 

The use of an iterative optimization method for obtaining the minimum of (16) requires a setup value 
from which to start iterating. An initial estimation of X£ close to its correct value reduces the probability 
to fall in a local extrema of the cost function. 

The method proposed in this paper to give an initial value for X£ consists of a non-iterative and exact 
method, which gives directly an accurate solution for Xq in absence of noise in odometry values. This 
method is based on the assumption that the robot moves in a plane and thus only the angle over its Z axis 
is needed. As explained below, this assumption allows to solve the problem with a rank deficient Linear 
Least Squares approximation, which is solved using the Singular Value Decomposition of the system's 
matrix and imposing a constraint that ensures a valid rotation. Its development is briefly introduced next. 

For a point M % of robot's model and at time k of the initialization trajectory, the image measurement 
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y\ results from the following projection: 




\IK c (R c (R£R 0 M 1 + T 0 + R,T£) + T c 



(18) 



where the matrix R k and vector T k represent the rotation and position of the robot in the floor plane 
given by the odometry supposing that X 0 = ^0 0 0^ . The rotation matrix R 0 and the offset T 0 
correspond with the initial pose X 0 

( 
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Xq, Xq, a 0 ) in form of rigid transformation: 
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(19) 



V 0 / 



where R 0 is a non-linear function of the orientation angle. The expression (18) depends non-linearly of 
vector Xq and so a new parametrization is proposed jointly with a transformation which removes the 
projective parameter A^. 

The point M l is replaced by the rotated M x = RqM\ removing thus the product between unknowns. 
The orientation angle in the pose is substituted by parameters a = cos(ao) and b = sin(oto), with 
the constraint a 2 + b 2 = 1. Using the new parameterization the expression (18) is transformed in 
the following: 




\K c (R c (r£m Xq + To + L t a 
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(21) 



where x^.' A y x 2 k ' A are the two first coordinates of Tj^. 

The new unknown vector, which correspond to the new parametrization of X$, is: 

T 



*=(xl x 2 a b (M Xo ) T ,- 
The expression (20) is decomposed in the following elements: 

(Ul\ 



l x 0 j 



(22) 



vi. 



vi 



K(R c {R£M Xq + To + L t a) + T c ) 



(23) 



where U k , V k y S\ are linear in terms of $ but not in terms of y k . 



L* Vk $ + b\ 



Ql 



L^ + b 



(24) 
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The projective scale is removed from the transformation by using vector product properties: 

AA 



a; 



fui\ 

vi 



fui\ 



(ui\ 

vi 



v 1 / W W \ S U 

where two lineally independent equations are obtained for $. 



0 



W 



(25) 



vlSl 



VI = 0 



(26) 



-u\Sl + U l k = 0 

Using all measurements inside Y L , a lineal system of equations is obtained in terms of $: 



/ {v\L l 
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(27) 



\ U K°S K + °U K J 



It is straightforward to show that system of (27) has a single-parameter family of solutions. If $ 0 is a 
possible solution, then $ 0 + is a solution for any ip 6 R, with A: 

/ /o\ T /o\ r \ 



A 



\ 



0 



0 

\T C>X 3 J 



(28) 



/ 



In fact, if A is normalized, it matches up with the eigen-vector associated to the zero eigenvalue of 
matrix A 7 A. 

Using the constraint that a 2 + b 2 = 1, and the singular value decomposition of matrix A, an exact 
solution of system (20) is obtained. 

Once $ is available, the solution Xq is obtained by inverting the parametrization: 

a 0 = tan" 1 (a, 6) M { = R^M l 0 X 0 a = (xl x 2 a 0 (M X ) T ■•• (M W ) T ) T (29) 

This method, although exact, is prone to error due to odometry noise and does not benefit from the 
Maximum Likelihood metric. However it is valid as a method to give an initial value for Xq before using 
the iterative approach. 

3.6. Degenerate Configurations 

The kind of trajectory performed by the robot during initialization has direct influence in the solution 
of Xq. There are three kinds of movements that yields degenerate solutions: 

• Straight motion: there is no information about the center of rotation of the robot and thus the pose 
has multiple solutions. 
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Rotational motion around robot axis: the following one-parameter (i.e., n) family of solutions 
gives identical measurements in the image plane: 

/ 0 \ 



M l n = nM l + (n - 1) 



0 



T 0n = nT 0 + (n - 1) 



(t c (i)\ 

Tc(2) 
\ 0 / 



with T c = R T C T C . 



X a k(n) = (T 0n a 0 Ml 



■■■ M. 



N 



(30) 



(31) 



Circular trajectory: under a purely circular trajectory the following one -parameter family of 
initialization vectors gives identical measurements in the image plane: 



M* = nM i + Rl T 0n = nTo + (n - l)R T c T c 



(32) 



X a k (n) = (T 0n a 0 Ml ■■■ M") 



(33) 



In practical cases it has been proved to be effective enough to combine straight trajectories with 
circular motion to avoid degeneracies in the solution of Xg . 

3. 7. Obtaining the Gaussian Equivalent of the Solution 

Once the minimum of (16) is reached we suppose that the resulting value of X$ is the mean of the 
distribution Xq. This section describes how to also obtain the covariance matrix. 

The covariance matrix Eg of the optimized parameters is easily obtained by using a local 
approximation of the term Yl — Yl in the vicinity of the minimum Xq using the following expression: 



^0 



(34) 



where J is the Jacobian matrix of Y with respect to parameters Xg . The Jacobian is available from 
the optimization method, in which is used to compute the iteration steps. 

4. Online Robot Localization 



In this section the solution to X% given the last pose information is derived. The fact that last frame 
information is available and the assumption of soft motion between frames allows to greatly simplify 
the problem. 

A special emphasis is given to the fact that any process handled by the system is considered a random 
entity, in fact a Gaussian distribution defined at each case by its mean vector and covariance matrix. The 
problem of obtaining pose and structure, encoded in X% given image observations and the previous 
time estimation X£_ x is viewed from the point of view of statistical inference, which means searching 
for the posterior probability distribution p(X%\Yi, ■ ■ ■ , Y^). That distribution gives the best estimation 
of XI given all the past knowledge available. In Figure 6, a brief overview of the online method 
is presented. 
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Figure 6. Overview of the online algorithm. 
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The online approach is divided into three steps: 

• Estimation Step: using the previous pose distribution p(X^_ 1 \Yi, ■ ■ ■ , Y k _i), defined by its mean 
X£_ x and covariance matrix and the motion model g a , a Gaussian distribution which infers 
the next state is given p(X%\Yi, ■ ■ ■ , Y^-i). 

• Robust Layer: the correspondence between image measurements and the 3D model of the robot 
easily fails, so a number of wrongly correspondences or outliers pollute the measurement vector 
Yfc. Using a robust algorithm and the information contained in the state vector, the outliers are 
discarded before the next step. 

• Correction Step: using an outlier-free measurement vector, we are confident to use all the 
information available to obtain the target posterior distribution p(X%\Yi, ■ ■ ■ , Yk) 

In all three steps it is required to propagate statistic processes over non-linear functions (/ and h). 
As was stated in the previous section we show how to face the problem using first order expansions as 
it offers more compactness and is more readable. As a consequence the "Estimation" and "Correction" 
steps are solved using the so called Extended Kalman Filter (EKE) equations, which have been already 
implemented in problems of similar complexity such as visual Simultaneous Localization and Mapping 
(SLAM) [5]. 

4.1. Estimation Step 

The estimation step uses the motion models available to infer the next pose of the 
robot which implies an increment in uncertainty. Starting from the last pose distribution 
p(X^_ 1 \Yx, ■ ■ ■ , Yk-x) = N(X^_ V the motion model g a and the noise included in odometry 

values, the following update is obtained: 

x k\k-i = 9 a { X kiUk) 

^fcifc-i = Jx^k-\Jx + Ju^wJu-, (35) 
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where X^ k _ 1 and S^ fe l are the mean and covariance matrix of distribution: 

P (xz\Y lr .. ,y fc _o 

The matrices Jx and J[/ are the first derivatives of the function g a with respect to X k _ x and Uk 
respectively. Usually J x in odometry systems is the identity, therefore, at this step, the covariance 
matrix S^fc-i resu l ts t0 be bigger in terms of eigenvalues, which means uncertainty. 

4.2. Correction Step 

The correction step removes the added uncertainty in the estimation by using image measurements. 
It passes from the distribution p(X%\Yi, ■ ■ ■ , Yfc_i) to the target distribution p(X%\Yi, ■ ■ ■ , Yk), which 
includes the last measurement. 

Using the estimation shown in (35), and knowing the correspondence between measurements with 
the camera and structure point of the state vector, the estimated measurement is given: 

Y^-i = h a (X a k ) (36) 
s y fc | fc _i = Jh^k\k-iJh + ^v (37) 
S X ay = £ Jk [ fc _ 1 «/ ft (38) 

where J h is the Jacobian matrix of the function h a with respect to X% and £y is block diagonal matrix 
with T. v on each block. Function h a performs the projection in the image plane of the camera of all 
visible points that form up the measurement vector The correction step itself is a linear correction of 
X^ k _ 1 and by means of the Kalman gain K G : 

K G = Ex-yE^ (39) 

X a k = -X"fe|fc-i + KciYk — Y k \k-i) (40) 

s fc = s fe|fc-i - K G T,xa Y (41) 

As it is stated in (41) the resulting is reduced compared to E^ fc l which means that after the 
correction step, the uncertainty is "smaller". 

4.3. Robust Layer 

The robust layer has the objective of removing bad measurements from Yk to avoid inconsistent 
updates of X% in the correction step. In this paper we propose to include the Random Sample Consensus 
(RANSAC) algorithm [22] between the estimation and correction step of the filter. The general idea is 
to found among the measured data Yk a set which agrees in the pose Xk, using the 3D model obtained 

inX%_ v 

The interest of applying RANSAC in a sequential update approach resides on several reasons: firstly 
it allows to efficiently discard outliers from Yk preventing algorithm's degeneracy, which happens even 
if the motion model is accurate. Secondly, compared to online robust approaches, where a robust cost 
function is optimized, RANSAC allows not to break the Kalman filter approach, as it only cleans the 
measurement vector of outliers. Furthermore we have observed experimentally that the RANSAC 
algorithm can be very fast between iterations, as only a few outliers are inside the data. (We use 
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the RANSAC implementation described in [17], which implements a dynamical computation of the 
outlier probability). 

The RANSAC method proposed in the commented framework obtains the consensus pose X k from 
the set of measurements Y k and the 3D model available in using the algorithm presented in [18]. 
For a robot moving in a plane, as it is the case with the mobile robot considered in this paper, it is 
enough to use a minimum of 2 correspondences between the model and the measurement which makes 
the RANSAC very fast for removing outliers. 

4.4. Image Measurements 

Contrary to the initialization case, in this step we have an accurate prediction of the tracked points 
in the image plane at each time instant, namely vector Using such prediction we can easily 

match the 3D points in the state vector with measurements taken in a measurement set T k , using the 
SIFT detector applied to current image I k . The matching is done in terms of distance in the image 
plane. Let y l k \ k _ 1 a feature inside vector Y^k-i and y{ a candidate obtained using SIFT method. We 
conclude that they are matched if \y J k — 2/L fe _ 1 |if < T max , where \ \m states for the Mahalanobis distance 
using the covariance T, yklk l computed from the matrix Sy fc|fe _ r The Mahalanobis distance allows to 
take into account the uncertainty predicted for y l k \ k _ 1 and also helps to select a threshold T max with a 
probabilistic criterion. 

5. Results 

This section describes the experimental setup developed to support the theoretical algorithms 
proposed in this paper. The experiments are divided in those performed over synthetic data and those 
run in a real implementation of the Intelligent Space in the University of Alcala (ISPACE-UAH) [23]. In 
both kind of experiments the same camera parameters are used, derived from the real device used in its 
real placement. The single camera consists of a CCD based sensor with resolution of 640 x 480 pixels 
and a physical size of 1/2 (around 8 mm diagonal). The optical system is chosen with a focal length of 
6.5 mm which gives around 45° of Field of View (FOV). The camera is connected to a processing and 
acquisition node through a IEEE 1394 port, which support 15 fps in RGB image format acquisition. The 
intrinsic parameters are the following: 

/„ = 636.7888 /„ = 637.5610 u 0 = 313.3236 v 0 = 210.6894 (42) 

The camera is placed with respect a global coordinate origin, as it is displayed in Figure 7. Camera 
calibration is performed previously to the localization task, using checkerboards as calibration patterns 
and the "Matlab Calibration Toolbox" implemented by [24]. The distortion parameters of the camera 
are not considered in this case. As it can be shown in Figure 7, the camera is placed in oblique angle, 
which is specially useful for covering large areas with less FOV requirements (less distortion) compared 
to overhead configurations. 
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Figure 7. Geometric distribution of the camera and robot's trajectory. 
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The robotic platform is connected to the same processing node controlling the camera by means 
of a wireless network. The camera acquisition and the odometry values obtained from the robot are 
synchronized. The control loop of the robot is internal, and it is prepared to follow position landmarks 
based on odometry feedback at faster sampling frequency than the image localization system (15 fps). 
Therefore, for each image acquisition cycle, the odometry system obtain several readings that are used 
by the internal loop control. In this paper the localization information provided by the cameras is not 
included in the control loop of the robot. 

5. 1. Simulated Data 

In this experiment the robot structure is simulated with a set of N = 10 points distributed randomly 
inside a cylinder with radius R = 0.5m and height h = lm. The odometry system is supposed to have an 
uncertainty a\ { = 10 and cr^ = 1 in linear and angular speed respectively. The initialization trajectory 
is shown in Figure 8. The image measurements are considered polluted with a Gaussian process of 
a"l = 10, independently of each measurement and image coordinate. 

Figure 8. Initialization trajectory of the robot in the experiment based on synthetic data. 
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To compare the performance of the initialization method proposed in Section 3., the following error 
magnitudes are described: 



cm 



£f =1 ||M*-M*J|2 



Ef=il|M^|| 2 



— | \Tq — To t gth\\ € a — I |ao — (X0,gth\ 



(43) 



where T 0)9 th, oco j9 th and M % th correspond to the ground truth values of the initialization vector Xg. The 
following two experiments are proposed: 

• Initialization errors in function of the odometry error. In this experiment the value of a\ ri and o\ 
are multiplied by the following multiplicative factor: 

p e (o.Ol 0.1 0.5 1 5 10) (44) 

The different errors of Equation 43 can be viewed in Figure 9a-c in terms of p. As it can be 
observed in the results, the M.C. (Complete correlated matrix H£) method outperforms the rest 
of approximations of S L , specially the full diagonal method M.I., which means that the statistical 
modelling proposed in this paper is effective. 

Figure 9. Experiment showing the different initialization errors in function of the amount of 
error in odometry readings. 
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Initialization errors in function of the trajectory length. The trajectory used in the experiment and 
displayed in Figure 8 is uniformly scaled by parameter p t , 



Pt e (0.2 0.4 0.6 0.8 1 1.2 1.4) 



(45) 



so that it can be guessed the relationship between trajectory length and initialization errors. In 
Figure 10 the initialization errors are displayed versus p t . 

In light of the results shown in Figure lOa-c, the M.C. method is capable of reducing the error no 
matter how large is the trajectory chosen. However, in the rest of the approximations of £ £ there is 
an optimal point where the initialization errors are minimum. This results make sense, as without 
statistical modelling large trajectories contain accumulative errors which usually affects the final 
solution of Xq. 
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Figure 10. Experiment showing the different initialization errors in function of the trajectory 
length performed by the robot. 



























I 


— M.C. 




M.B. 










-♦-M.D 
—M.I. 



















Value of factor p t for scaling the trajectory 

(a) e M vs p t 





-e-M.C. 








M.B. 








-♦-M.D 








—M.I. 











































—M.C. 












M.B. 
-♦-M.D 












—M.I. 



































Value of factor p } for scaling the trajectory 

(b) e T vs p t 



Value of factor p ( for scaling the trajectory 

(c) e a vs p t 



Both experiments clearly manifest that the complete matrix £ L , with all its cross-correlated terms 
(M.C), outperforms the rest of proposals, especially when it is compared to the case where £l is the 
identity matrix (M.I.), which means no statistical modeling. 

5.2. Real Data 

The initialization experiment proposed in this paper consists of a robot performing a short trajectory 
from which its 3D model and initial position is obtained using the results of Section 3.. We present a 
comparison of the initialization results using three different trajectories. Each one of the trajectories 
is displayed in Figure 11 as a colored interval of the whole trajectory performed by the robot in 
the experiment. 

Figure 11. Intervals of robot's trajectory used for initialization. 
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The results of the initialization method on each of the 3 trajectories selected is shown in Figure 12. 
Depending on the trajectory used, the features viewed by the camera vary and thus the corresponding 
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initialized 3D model. As it can be seen in Figure 12, the 3D model is accurate and its projection matches 
with points in the robot's structure in all cases. It must be remarked that on each case, the position 
obtained as a result of the initialization {i.e., X 0 ) corresponds to the first position of each interval. 

Figure 12. Single camera initialization results. On each row it is shown the resulting 
reconstruction and its projection in the image plane of the camera. 




The sequential algorithm is tested using the whole trajectory shown in Figure 1 1 . The initial pose and 
3D model are the result of the initialization results shown in Figure 12e,f. The estimated position of the 
robot, compared to a "ground truth" measurement, is presented in Figure 13a,b. 
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Figure 13. Online robot's pose compared with the "ground-truth". 




The ground truth is obtained by means of a manually measured 3D model of the robot. The model is 
composed of 3D points, that are easily recognized by a human observer. By manually clicking points of 
the 3D model in the image plane, and by using the method proposed in [18], the reference pose of the 
robot is obtained in some of the frames of the experiment. 

Another experiment is proposed to test the online algorithm with occlusions and a larger path followed 
by the robot. In Figure 14a it is shown the geometric placement of the camera and in Figure 14b it is 
shown the geometric model obtained during initialization. 

Figure 14. Scene geometry and 3D model obtained during initialization. 




(a) (b) 



The resulting localization error is shown in Figure 15b. In Figure 16 it is shown several frames, 
indexed by the time sample number k, presenting hard occlusions between the camera and the robot 
without losing the tracking. The RANSAC method used in the Kalman loop avoid erroneous matches in 
the occluded parts of the object. Besides, in those frames with completely occluded features, only the 
estimation step of the Kalman filter is done, which is accurate enough for short periods of time. 
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In both sequential experiments the 3D model is composed of 8 to 10 points and no extra geometrical 
information is introduced in the state vector. As a future line of study a simultaneous reconstruction and 
localization approach can be adopted to introduce online extra 3D points in the state vector as the robot 
is tracked. A very similar approach is done in the Simultaneous Localization and Mapping (SLAM) 
problem, and some of their solutions [5] are perfectly applicable to the problem assessed in this paper. 
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6. Conclusions 

This paper has presented a complete localization system for mobile robots inside intelligent 
environments with a single external camera. The objectives of obtaining the pose of the robot based 
on its natural appearance has been tackled in the paper using a reconstruction approach followed by a 
sequential localization approach. 

The main contributions of this paper are summarized in the following list: 

• The initialization step provides a non-supervised method for obtaining the initial pose and structure 
of the robot previously to its sequential localization. A Maximum Likelihood cost function 
is proposed, which obtains the pose and geometry given a trajectory performed by the robot. 
The proposal of this paper allows to compensate for the odometry drift in the solution. Also 
an exact initialization method has been proposed and the degenerate configurations have been 
identified theoretically. 

• The online approach of the algorithm obtains the robot's pose by using a sequential Bayesian 
inference approach. A robust step, based on the RANSAC algorithm, is proposed to clean the 
measurement vector out of outliers. 

• The results show that the proposed method is suitable to be used in real environments. 
The accuracy and non-drifting nature of the pose estimation have been also evaluated in a 
real environment. 

The future research must be oriented to scale the problem for large configurations of non-overlapped 
cameras and multiple robots, where extra problems arise, such as to automatically detect the transitions 
between cameras and online refinement of the geometric models. It is important from our point of view 
to, in the future, combine the information given by the system proposed in this paper with information 
sensed onboard the robots using cameras. This approach allows to jointly build large maps attached 
to information given by the cameras, so that robots can be localized and controlled to perform a 
complex task. 
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